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Using D-theory we construct a new efficient cluster algorithm for the Ising model. The con- 
struction is very different from the standard Swendsen-Wang algorithm and related to worm al- 
gorithms. With the new algorithm we have measured the correlation function with high precision 
over a surprisingly large number of orders of magnitude. 
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1. Reformulation of the Ising Model 

D-theory [1] is an alternative way of regularising field theories by dimensional reduction of 
systems involving discrete variables, that allows us to develop cluster algorithms. One can also use 
this method to describe and simulate a simple spin system such as the Ising model. Besides being 
interesting in its own right, we hope to gain some insights that may be useful for the construction 
of cluster algorithms for gauge theories. 

We consider the Ising model with the classical Hamilton function J^[s] = —JT,x,fi^x^x+fi^ 
and rewrite it as a quantum system with the Hamiltonian H = —JY,x.ti ^x^x+fi- ^'^ trace- 
structure of the partition function, Z = Trexp(— j8//), we can perform a unitary transformation 
rotating to a basis containing a^: 
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We then perform a Trotter decomposition of the Hamiltonian. In one dimension this implies the 
following splitting into two parts. 



X odd 



while in two dimensions one needs four parts 
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The partition function can then be turned into a path integral with M time steps of width £ (j8 = 
eM). The Euclidean time serves as an additional dimension whose extent j8 determines the inverse 
temperature. For example in the 1 -dimensional Ising model, the Trotter decomposition of the 
partition function takes the form 



Z = Trexp(-e//i - sHx)" = Trexp(-£7?i)*' exp{-eH2)'^ . 
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In this way we obtain a checkerboard decomposition where the action resides on the shaded pla- 
quettes (see figure [ij). 
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Figure 1: Trotter decomposition in one dimension with a spin-configuration: crosses and cuxles represent 
the two spin states. 

In general the Trotter decomposition is affected by an order error. However, since in the 
Ising model H\ and H2 commute, it is an exact rewriting of the partition function and even with the 
minimal number of time-slices one obtains the continuous time limit. 
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2. Construction of the Algorithm 

In our basis the transfer matrix of a single plaquette T = e.xp{—eG^Gl^^) takes the form 
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with c = cosh(£/) and s = sinh(£/). From this we can see that there are only eight physically 
allowed plaquette configurations (see figure ^. 

In order to imply constraints on the relative orientation of the spins on a plaquette, we propose 
plaquette breakups that bind spins together. These constraints must be maintained in the cluster 
update. The breakups give rise to clusters and these clusters can only be flipped as a whole. In this 
way, by construction, one indeed maintains the constraints implied by the cluster breakups. We 
choose to use the A, Bi, and breakups shown in figure |3[ 
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Figure 2: Allowed plaquettes. 



Figure 3: Cluster breakups. 



These breakups decompose the transfer matrix as follows: 



T =A 




+ B 




+ B2 




(2.2) 



From eq.([2. 1|), one obtains Bi = B2 = s and A = c — s. Using these cluster breakups we can now 
simulate the system. In figure ^ we show an example of a cluster update in which the dashed cluster 
has been flipped. 
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Figure 4: A configuration with A, By, and B2 breakups with a multi-cluster update. 

The energy density is related to the fraction P of interaction plaquettes with a transition (spins 
opposite in time on a plaquette) by the equation 



{H) = -J{t+^^—^{P)), 



(2.3) 
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where t = tanh(£7). 

Analysing the average cluster-size of the algorithm in multi-cluster mode shows that it is in- 
deed very different from the Swendsen-Wang algorithm [2]. In two dimensions, with the Trotter 



decomposition of eq.(|1.3[) and the minimal number of four time-slices, we observe a peak at the 
critical temperature, whereas Swendsen-Wang clusters still grow in the broken phase. This is shown 
in figure 0. 



Figure 5: Average cluster-size of the new algorithm (diamonds) compared to the Swendsen-Wang algorithm 
(crosses) on a 100^ lattice. 



3. Measurement of the 2-Point Correlation Function 

The correlation function is defined in our basis as (ajoj). The matrices are off-diagonal 
and flip the spins they act on. Thus they can be viewed as violations of spin conservation, i.e. two 
opposite spins reside on the same site (see figure 
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1 

Figure 6: The action of the operator a' on a site, seen as a violation of spin conservation. 

What we do is the following: we introduce two violations on a random site (where they cancel 
each other), build the clusters, and flip them with probability p = ^, unless it is a cluster with one 
or two violations. There we flip a fraction of the cluster in order to move the violations to different 
positions. For example in figure^ the dashed part of a cluster with violations has been flipped. 

When the violations are on the same site, we introduce them randomly somewhere else. By 
histogramming the spatial distance between the two violations, we get a very accurate measurement 
of the correlation function. One then normalises the correlation function to 1 at zero distance. In 
figure^ this is compared to the analytic result [3]. 
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Figure 7: Partially flipping a cluster with violations (circles). 

4. Improved Estimator for the Susceptibility 

The susceptibility can be measured by summing over the correlation function. This is already 
an improved estimator because one gets only positive contributions. One can also measure the 
average frequency of reaching spatial distance zero. This corresponds to the inverse susceptibility. 
In a large system at high temperature this is not too efficient yet and can be further improved. 

We project the clusters with violations to the spatial volume and measure their overlap. This 
gives us the probability that two violations are at zero spatial distance in the next cluster update 
and, again, corresponds to measuring the inverse susceptibility. In this way we gain statistics by a 
factor of the cluster-size squared |Cp with an effort proportional to the cluster-size |C|. 



5. Worm-Algorithm 

Instead of using clusters, we can — with the same breakups — move the two violations locally 
with a worm-type Metropolis algorithm [4]. 
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Figure 8: Worm moving the violations and thus updating a given configuration. 

Interestingly we can also measure n-point functions by introducing n violations. These viola- 
tions are moved in the same way as above. 

5.1 Snake- Algorithm — Improved Correlation Function Measurement 

An efficient trick [5] allows us to measure an exponentially suppressed signal with a linear ef- 
fort. The correlation function is exponentially suppressed in the high-temperature phase and is thus 
hard to measure at large distances. The correlation function is the ratio of two partition functions 
and it can be rewritten as the product of many fractions which can be measured individually: 

Z(x) Z(1)Z(2) Z(x) 



z(o) z(o)z(i) ■■■z(x-i)' 



(5.1) 
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In figure [T^ we compare the numerical results at very high temperature with the analytical expres- 
sion [3]. 




Figure 9: Correlation function on a 64^-lattice at Figure 10: Correlation function on an 80^-lattice at 
j3 = 0.42 (crosses) and j3 = 0.46 (diamonds) com- j3 = 0.01. 
pared to the analytic result. 



6. Conclusion 

Using D-theory we have constructed a new efficient algorithm for the Ising model. We have 
an improved estimator for the susceptibility. Combining the algorithm with the snake algorithm 
one can measure the correlation function in regions otherwise hardly accessible. We would like to 
apply this approach to other theories, possibly gauge theories. 
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